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1  IMTRQDUCTION 

The  wave  drag  of  slender  pointed  bodies  of  general  cross-section  has  been  shown^t 
within  the  restrictive  assumptions  of  slender-body  theory,  to  be  given  by 


D  •  -  4itp 


f'(x)f'(s)  log  lx  -  s|  dxds 


+  4irof(il)  I  f'(x)  log(e  -  x)  dx 


2  2 
1 6ir  pU 


f  If 


do  0(t  log  t) 


where  x  is  the  streamwise  axis,  f  is  the  strength  of  a  smooth  source  distribution 
along  the  x  axis  (for  0  <  x  <  i) ,  I.  is  the  body  length,  t  is  a  thickness  paramc  -.er 
and  p  and  U  are  the  free-stream  density  and  velocity  respectively.  The  third  term 

in  equation  (1)  is  a  contour  integral  round  the  base  of  the  body,  vdiere  i)>  is  the  slender- 

body  potential  and  v  is  the  normal  to  the  contour  in  the  base  plane. 

The  source  distribution  f  can  be  related  to  the  cross-sectional  area  of  the  body 
by  the  simple  relation 

f(x)  =  --^SMx)  (2) 

permitting  the  wave  drag  of  equation  (I)  to  be  expressed  directly  in  terms  of  the  geometry. 

For  non-lifting  configurations  which  have  a  pointed  base  or  which  merge  smoothly  into 

a  cylinder  parallel  to  the  free  stream,  the  second  and  third  terms  on  the  right-hand  side 
of  equation  (I)  are  zero,  and  equations  (I)  and  (2)  reduce  to  the  more  familiar  Von  Karman 
wave  drag  integral  form: 


-i/  /=■• 


(x)S"(8)  log|x  -  s|  dxds  . 


The  first  practical  application  of  equation  (3)  to  non-slender  bodies  was  made  by 
,  2 

Whitcomb  ,  who  observed  that,  at  transonic  speeds,  the  drag  of  a  wing-body  combination 
was  approximately  the  same  as  that  of  an  equivalent  body  of  revolution  with  the  same 
cross-sectional  area  distribution.  Although  Whitcomb  stated  his  area  rule  in  terms  of 
areas,  for  theoretical  justification  it  is  necessary  to  consider  the  body  in  terms  of 
source  distributions.  The  equivalent  axial  body  is  then  one  which  has  a  source  distribu¬ 
tion  formed  from  the  sum  of  the  sources  lying  in  planes  perpendicular  to  the  free-stream 
direction  (which  are  approximately  Mach  planes  in  transonic  flow) . 


When  the  Mach  number  is  not  close  to  1,  the  Mach  lines  are  not  nearly  perpendicular 
to  the  stream  and  some  modification  to  the  area  rule  is  required.  In  supersonic  flow 
the  equivalent  source  strength  is  taken  as  the  sum  of  the  sources  lying  on  the  Mach 
plane 

X  -  By  sin  8  -  Bz  cos  8  -  Xq  (4) 

where  8  (the  cylindrical  polar  angle)  and  Xq  (the  intersection  of  the  plane  with  the 
streamwise  body  axis),  are  parameters  \diich  define  the  particular  Mach  plane.  The  wave 

3 

drag  of  a  non-lifting  body  can  then  be  expressed  as 

2ir  1(6)  1(6) 

f (x,6)f ' (s,6)  log|x  -  s|  dxdsd6  .  (5) 

The  evaluation  of  this  triple  integral  can  be  simplified  if  the  wing  sources  (f^) 
and  the  fuselage  sources  (f^)  are  separated.  That  is,  writing  f'(x,6)  as  f^  +  f^  . 
equation  (5)  splits  into  terms  representing  the  drag  of  the  fuselage  alone,  the  drag  of 
the  wing  alone  and  the  wing-fuselage  interference  drag.  If  the  fuselage  is  slender,  the 
drag  of  the  fuselage  can  be  obtained  from  equation  (3)  using  normal  cross-sectional 
areas  (which  do  not  vary  with  8),  and  also  the  order  of  the  integration  of  the  inter¬ 
ference  term  can  be  reversed  allowing  the  8  integration  to  be  simplified.  If  the  wing 
is  described  adequately  by  sources  lying  in  a  plane,  the  expression  for  the  drag  of  the 
wing  alone  can  also  be  simplified  analytically,  although  for  the  wing  in  isolation  a 
better  estimate  of  the  drag  may  often  be  obtained  by  empirical  methods. 

When  applied  to  aircraft  which  really  are  slender  (such  as  Concorde)  the  above 
.  4 

techniques  give  good  estimates  of  the  wave  drag  .  More  recently,  however,  various 
attempts  to  estimate  the  wave  drag  of  fighter-type  aircraft  have  had  mixed  success,  with 
sometimes  significant  differences  occurring  in  the  estimates  from  methods  which  only 
appear  to  differ  in  the  details  of  their  application.  It  is  to  be  expected  that  the 
application  of  the  methoas  to  aircraft  which  are  not  slender  would  lead  to  decreased 
accuracy;  fdiat  has  proved  to  be  disappointing,  however,  is  the  lack  of  pattern  in  the 
error  and  the  consequent  limited  success  of  modifications  intended  to  reduce  this  error. 

One  potential  source  of  difficulty  for  some  methods  is  the  treatment  of  the  fuse¬ 
lage  and  wing  by  differing  techniques  (namely  normal  and  oblique  cuts  respectively) . 

Modem  fighter  aircraft  often  have  the  root  of  the  wing  extended  forward  along  the 
fuselage  or  blended  into  the  fuselage,  so  that  it  is  not  clear  how  the  components 
should  be  distinguished.  A  conceptually  simpler  technique,  which  removes  this  difficulty, 
is  to  generate  the  area  distributions  for  the  triple  integral  using  'all  oblique  cuts', 
that  is,  to  obtain  the  cross-sectional  areas  from  the  intersection  of  the  wing  and  fuse¬ 
lage  with  the  Mach  planes  given  by  equation  (4) ,  and  to  evaluate  the  triple  integral 
directly. 

Which  of  the  two  approaches  is  the  more  applicable  to  less  slender  bodies  is  unclear, 
but  the  latter  approach  would  appear  to  be  at  least  worthy  of  investigation.  The  reason 


for  its  neglect  has  been  mechanical  rather  than  theoretical,  in  that  it  has  proved  to  be 
surprisingly  difficult  to  find  the  required  oblique  cross-sectional  areas  of  even  quite 
simple  bodies. 

In  section  2  of  the  Report  a  novel  numerical  technique  is  presented  for  obtaining 
oblique  cross-sectional  areas  for  general  body  shapes  and  in  section  3  the  evaluation 
of  wave  drag  from  the  lengthwise  distribution  of  these  areas  using  Eminton's^  method  is 
briefly  described.  A  computer  program  listing  with  subroutines  for  several  classes  of 
bodies  is  given  in  the  Appendix  (examples  from  these  classes  are  used  in  section  A), 
together  with  a  description  of  the  input  parameters.  The  accuracy  of  prediction  of  wave 
drag  for  various  bodies  is  discussed  in  section  4. 

2  BODY  DESCRIPTION 

The  numerical  description  of  the  body  needs  to  be  both  efficient  in  its  use  of 
computer  storage  and  of  such  a  form  that  oblique  cross-section  areas  may  easily  be  found. 
In  this  section  a  body  description  meeting  both  these  requirements  is  developed. 

Consider  a  line  parallel  to  a  convenient  body  axis,  x  ,  which  intersects  the  body 
as  shown  in  Fig  1 .  The  line  is  uniquely  identified  by  its  y  and  z  coordinates  (y^ 
and  Zq  say)  and  the  section  of  the  line  inside  the  body  is  defined  by  the  x  values 
(x. ,  with  X.  ,  >x.)  at  its  intersections  with  the  body.  For  most  simple  bodies  two 
X  values  are  sufficient  (X|  and  x^),  but  for  more  complex  bodies  multiple  intersections 
may  occur  requiring  further  x^^  values. 

Using  a  slender-body  analogy,  the  body  volume  local  to  the  line  y^,  z^  is  con¬ 
sidered  to  be  concentrated  on  that  line.  The  area  of  an  oblique  cross-section  of  this 
element  can  then  immediately  be  related  to  the  normal  cross-section  area  as  in  slender- 
body  theory.  Thus  if  A(x)  is  the  elemental  cross-sectional  area  associated  with  y^, 
(corresponding  to  S  used  in  equation  (2)  for  slender  body  theory),  then  close  to  the 
intersection  points  x^  ,  A  will  vary  according  to  the  slope  of  the  body,  and  remote 

from  these  points  it  will  be  constant.  If  the  constant  cross-section  area  is  A 

max 

inside  the  body  then  the  element  can  be  adequately  defined  by  y^,  Zg,  x^,  ®nd  the 

variation  of  A  near  the  intersection  points  x.  . 

To  achieve  a  compact  representation  of  the  body  in  the  Computer  it  is  necessary  to 
reduce  to  a  minimum  the  number  of  values  defining  each  element.  This  is  achieved  here 
by  considering  the  computer  representation  and  the  physical  representation  of  the  data 
together.  In  the  computer  it  is  natural  to  locate  the  numbers  defining  the  elements  in 
rectangular  arrays.  The  y^  and  Zg  arrays  can  be  discarded  if  the  physical  position 
of  the  elements  in  space  can  be  simply  related  to  the  position  of  numbers  in  the  arrays 
which  contain  x^  and  the  other  numbers  necessary  to  describe  the  element.  The  simplest 
mapping  in  which  the  y^,  Zg  distribution  is  a  rectangular  array  of  points  in  physical 
space  is  wasteful,  because  some  of  the  (yg,Zg)  lines  may  not  intersect  the  body.  In  Fig  2 
showing  a  front  view  of  the  body,  it  is  indicated  how  the  physical  representation  of  a 
typical  column  of  Zg  points  is  compressed  such  that  nearly  all  the  lines  (Yqi^q)  inter¬ 
sect  the  body.  The  extreme  points  (denoted  by  I^  -  I  and  I^  »  N^)  are  arranged  to  lie 
within  Az/2  of  the  surface,  where  Az  is  the  distance  between  adjacent  points  and  is 
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independent  of  z  .  Note  that  and  represent  the  outer  limits  of  the 

elemental  area  representation  (rather  than  the  z  coordinates  of  the  end  points  of 
the  column)  such  that 

z  -  z  .  “  N  Az 

max  min  z 


(6) 


and 


z  .  +  (I  -  0,5)Az 

min  z 


I 

z 


1.2 . .  (7) 


In  this  mapping  the  total  number  of  points  in  each  column  is  kept  constant  (at  N^)  but 

the  z  .  .  z  and  Az  values  vary  from  column  to  column.  The  y  spacing  of  the 

min  ’  max 

columns  of  y^,  z^  values  is  assumed  to  be  constant,  so  that  the  interval  between  the 
columns.  Ay  ,  is  constant. 

The  body  is  also  assumed  to  be  symmetric  about  the  plane  y  =  0  such  that 


and 


yo^V  ° 


I 

y 


1,2,...,N 


(8) 


Ay 


y  /N 
■'max  y 


(9) 


With  this  representation,  the  value  of  is  given  by  AyAz  and  the  definition  of  an 

element  needs  only  its  x^^  values  and  information  on  the  variation  of  A  near  the  inter¬ 
section  points  (x  =  X. ) .  The  variation  of  A  near  the  intersection  points  depends  on 
the  surface  slope,  which  can  be  approximated  from  the  values  of  x.  for  the  adjacent 
elements.  Because  the  physical  array  of  (yQ,Zp)  points  is  not  rectangular,  the  only 
immediately  useful  adjacent  x^  values  in  the  array  are  those  in  the  same  column 
(ie  in  the  adjacent  z  positions);  the  x^  values  in  adjacent  y  positions  do  not 
correspond  to  the  same  value  of  z  and  use  of  them  would  require  interpolation.  There 
is  an  additional  requirement  on  the  variation  of  A  with  x  ,  in  that  the  elemental 
bodies  are  assumed  to  be  slender,  requiring  an  upper  limit  on  the  permitted  value  of 
dA/dx  .  With  these  restrictions,  the  method  used  to  approximate  dA/dx  is  shown  dia- 
gramatically  in  Fig  3.  The  variation  in  the  area  distribution  A  between  AyAx  and 
zero  is  taken  to  be  the  normal  cross-sectional  area  of  a  tapered  wedge  with  a  leading 
edge  lying  in  a  vertical  plane  y  -  constant  ,  and  oblique  to  the  x-axis.  The  wedge  is 
tapered  over  a  distance  of  2Ay  to  satisfy  the  slenderness  requirement  and  consequently 
the  maximum  value  of  dA/dx  is  Az/2  .  The  angle  of  the  slant  of  the  leading  edge  of  the 
wedge  is  obtained  from  the  x.  values  of  the  adjacent  Zg  points.  This  approximation 
to  dA/dx  is  very  easy  to  calculate  and  in  practice  proves  to  be  an  adequate  representa¬ 
tion  of  the  change  in  the  area.  The  representation  of  each  element  can  now  be  described 
by  its  intersection  points  only,  that  is  for  simple  bodies  the  entry  and  exit  points 
X|  and  X2  . 

In  the  computer  these  x  values  are  conveniently  stored  in  a  three-dimensional 

array,  which  in  the  program  as  listed  in  the  Appendix  is  called  the  XBM  array.  The  size 

of  this  array  is  N  ,  N  ,  IB  where  IB  is  the  maximum  number  of  intersections  (ie  the 
y  z 
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in 
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maximum  value  of  i).  The  computation  of  oblique  cross-sectional  areas  with  this  rep¬ 
resentation  is  now  no  more  complicated  than  that  of  normal  cross-sectional  areas,  for  it 
is  necessary  only  to  sum  the  contribution  from  each  element  at  the  x  value  where  its 
centreline  is  cut  by  the  plane  of  equation  (4).  To  find  an  oblique  cross-sectional 
area  where  the  normal  to  the  cutting  plane  makes  an  angle  with  the  element  axis,  the 
contribution  from  each  element  is  either  zero  or  4yAz  sec  i|)  ,  except  for  the  elements 
cut  near  their  end  points  when  the  correction  due  to  taper  is  applied.  This  simple 
repetitive  process  is  ideally  suited  to  solution  on  a  computer. 

3  EVALUATION  OF  THE  WAVE  DRAG 

Given  the  ability  to  calculate  oblique  cross-section  areas  at  specific  x  values, 
the  triple-integral  of  equation  (5)  can  be  evaluated  to  give  the  wave  drag  by  using 
Eminton's  method^.  This  method  depends  on  representing  the  area  distribution  as  a  sum¬ 
mation  of  sine  waves  (similar  to  a  Fourier  series),  which  on  substitution  in  the  integrand 
of  the  triple-integral  gives  terms  which  can  be  integrated  analytically  with  respect  to  x  . 
The  method  has  been  widely  used  and  is  not  discussed  further  here. 

The  computer  evaluation  is  performed  using  the  Fortran  subroutine  WAVEDR  listed 
with  the  main  program  in  the  Appendix.  As  this  subroutine  may  usefully  be  employed  other 
than  in  the  present  program  it  is  written  to  be  independent  of  the  main  program.  Its 
only  external  reference  is  to  the  subroutine  INVERT,  which  is  a  subroutine  to  invert 
a  matrix  A  onto  itself.  The  values  required  in  the  evaluation  are  passed  through  to 
parameter  list.  The  cross-sectional  areas  to  be  used  in  the  calculation  are  contained  in 
the  vector  AREA(NX) .  The  x  spacing  of  the  areas  (Ax)  is  given  by  DX  in  the  program, 
and  on  exit  from  the  subroutine  the  parameter  WDRAG  has  the  value  D(6)/q  in  the  units 
of  AREA.  In  the  main  program,  WAVEDR  is  used  to  obtain  D/q  at  NTHETA  values  of  6 
between  (and  including)  6-0  and  6  =  ir  .  The  body  wave  drag  is  then  taken  as  the 
appropriate  mean  of  these  values,  that  is  with  contribution  from  the  end  values  at  6  =  0 
and  9  “  TT  halved  and  the  mean  divided  by  ir  . 

4  ASSESSMENT  OF  THE  METHOD 

Wave  drag  estimates  from  the  supersonic  area  rule  method  discussed  here  are 
compared  with  more  accurate  wave  drag  estimates  to  validate  the  computer  program  and  to 
assess  the  error  involved  in  using  supersonic  area  rule  to  estimate  the  drag  of  non- 
s lender  bodies. 

4. I  Comparison  with  slender  body  theory 

To  investigate  the  errors  introduced  by  the  approximation  to  the  body  geometry 
and  the  other  approximations  used  in  the  present  method,  the  wave  drag  of  a  range  of 
Sears-Haack" ’  bodies  and  some  non-axi symmetric  bodies  is  computed  and  compared  with 
the  values  obtained  analytically  using  slender-body  theory. 

(a)  Sears-Haack  bodies 

The  shape  of  a  Sears-Haack  body  of  unit  length  is  given  by^ 


-  itR" 


128V 

Sir 


x^(l  - 


X)’ 


S(x) 


(10) 


The  wave  drag  of  thia  body  according  to  slender  body  theory  is  D 


where 


D  I28V^ 
q  “  K 


(11) 


Eliminating  the  volume  V  between  equations  (10)  and  (II)  we  obtain  an  equation  for  the 
body  shape  «rith  the  wave  drag  as  a  parameter,  ie 


The  substitution  of  a  value  of  D/q  in  this  equation  defines  the  geometry  of  a 
test  body  with  this  theoretical  value  of  D/q  .  Thus  for  this  body,  the  computer  program 
(with  M  ■  1)  ought  to  recover  the  prescribed  value  of  D/q  .  The  variation  in  the 


accuracy  of  the  program  with  the  thickness  of  the  body  and  the  number  of  elements  (N  ,N  ) 

y  ^ 

defining  it  is  shown  in  Figs  4  and  5  respectively.  In  Fig  4  the  error  obtained  when 


estimating  the  wave  drag  of  a  Sears-Haack  body  represented  by  2500  elements 


(.ie  N  »  N  =  50)  is  shown.  For  a  range  of  D/q  from  zero  to  0.1  (representing  a  body 

y  2 

thickness  to  length  ratio  from  0  to  0.327)  we  can  see  that  the  numerically  derived  value 


is  within  +1%  of  the  theoretical  value.  The  large  range  of  body  thickness  represented 
can  be  appreciated  from  the  body  side  views  shown  for  D/q  equal  to  0.02  and  0.08.  The 
pattern  of  fluctuation  in  the  error  is  of  little  significance  because  it  depends  on  the 
particular  relationship  between  the  elements  and  the  cross-sectional  cuts.  As  the 
number  of  elements  is  increased  it  is  shown  in  Fig  5  how  the  accuracy  increases,  as  might 


be  expected,  like  I/CN^N^). 


(b)  Norman's  'single  bulge'  body 


A  non-axis}niii>etric  body  for  which  oblique  cross-section  areas  have  been  derived 

8 

analytically  is  illustrated  in  Fig  6.  The  body  is  symmetrical  about  the  plane  y  =  0  and 
is  described  by  Norman  as  a  'single  bulge'  body.  The  body  cross-section  for  x  ^  0 
(ie  upstream)  is  square.  Between  x  =  0  and  x  =  6  the  body  cross-section  changes  smoothly 
to  an  'inverted  T'  section  in  such  a  manner  that  the  normal  cross-section  area  of  the 
body  remains  constant.  Downstream  of  this  change  (ie  x  >  6)  the  cross-section  shape  is 
again  constant.  The  body  is  defined  mathematically  by 


Top  surface 

z 

-  1.5 

(13) 

Bottom  surface 

z 

-  -  1.5 

(14) 

Sides  (z  >  0) 

y 

-  ±(I.5-H(x)) 

(15) 

Sides  (z  <  0) 

y 

-  ±(1.5+H(x)) 

(16) 

with  sufficient  of  the  plane  z  -  0  to  complete  the  surface.  The  function  H(x)  is 
defined  by 


H(x)  -  0 

H(x)  ■  i  -  J  cos(irx/6) 


X  <  0 
0  <  X  <  6 


(17) 

(18) 


This  body  critically  tests  the  accuracy  of  the  present  method  because  of  the  large 
areas  of  high  slope  {is  dy/dx)  on  its  vertical  slab  sides.  As  the  cross-sectional 
area  of  the  body  is  constant,  the  wave  drag  using  normal  cross-sectional  areas  in 
equation  (3)  is  zero.  The  wave  drag  from  oblique  areas  will  not  be  zero,  however. 


The  actual  values  of  the  drag  coefficient  found  from  the  analytically  derived  oblique 
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areas  are  shown  in  Fig  7. 


2  4 

By  introducing  the  supersonic-hypersonic  similarity  parameter  t(M  -  I)  ,  Fig  7 


can  be  made  to  represent  the  wave  drag  of  the  class  of  bodies  formed  by  a  linear 
stetching  of  the  body  defined  by  Korman  along  th^  x  axis.  The  function  H(x)  is 
then  defined  by 

H(x)  *  0 


H(x)  =  4  “  4  cos(Trx/6T) 

H(x)  =  1 


X  0 

0  <  X  <  6/t 
X  >  6/t 


(20) 

(21) 

(22) 


where  unit  value  of  t  recovers  Norman's  body  as  described  by  equations ( 1 3)  to  (19). 


The  difference  in  the  wave  drag  between  the  present  method  and  that  shown  in  Fig  7, 
is  indicated  in  Fig  8  for  M  =  1 .4  and  M  =  2.0.  To  obtain  a  1%  accuracy  it  is  necessary 
to  use  much  greater  detail  in  the  body  specification  than  for  the  Sears-Haark  body. 


values  of  200  being  required  for  and  .  The  reasons  for  this  are  twofold.  First 


the  body  has  slab  sides  with  large  gradients  near  x  =  3/t  which  require  detailed 
representation.  The  rod  geometry  described  in  Fig  3  for  the  area  distribution  at  the 
ends  of  the  rod  is  only  partially  effective  for  vertical  slab  sides.  This  makes  the 
geometry  of  this  particular  body  a  good  test  of  the  method  and  accounts  for  the  higher 


values  of  N 


The  second  reason  is  that  the  body  has  a  constant  cross-sectional  area 


normal  to  the  flow  direction  and  hence  the  wave  drag  is  due  to  interference  effects 


which  only  appear  when  M  >  1 .  This  interference  drag  is  much  smaller  than  the  drag 

g 

normally  associated  with  the  thickness,  for  example,  the  drag  coefficient  at  M  =  1  of 


the  top  and  bottom  halves  of  the  body  taken  separately  is  0,138,  which  is  more  than 
three  times  the  maximum  value  shown  in  Fig  7.  Thus  a  better  representation  of  the  body 
shape  than  for  a  closed  axisymmetric  body  is  required  to  obtain  the  1%  accuracy. 

4.2  Comparison  with  wave  drag  from  other  methods 

To  determine  whether  the  use  of  oblique  cuts  has  any  advantages  in  accuracy  over 
other  linear  methods,  it  is  necessary  to  estimate  the  wave  drag  of  bodies  whose  wave 
drag  is  known  accurately. 

(a)  Cone-cylinder  bodies 

Consider  axisymmetric  cone-cylinder  bodies  where  the  cylinder  extends  downstream 
to  infinity.  The  wave  drag  of  a  cone-cylinder  with  the  cross-section  area  of  the 


cylinder  as  reference  area,  is  the  same  as  the  constant  pressure  coefficient  on  the 


cone^.  In  Fig  9  the  wave  drag  of  cone-cylinders  with  10°  and  15°  semi -vertex  angles 


is  shown  for  a  range  of  Mach  numbers,  obtained  from  exact  theory  and  several  approximate 
theories*^.  The  exact  result  is  given  by  the  continuous  line.  The  results  from  slender 
body  theory  are  shown  by  the  fine  dotted  line  and  those  for  linear  theory  using  the  exact 
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boundary  conditions  on  the  cone  by  the  dot-dash  line.  Both  of  these  theories  are  clearly 
inadequate  for  this  case.  A  second-order  approximation  due  to  Broderick*'  and  shown  by 
the  plain  dashed  line  is  rather  better,  but  even  that  gives  an  error  of  over  lOZ  for  the 
15°  cone  above  M  -  1.5. 

The  Supersonic  Area  Rule  (SAR)  method  described  here  gives  the  'cross'  line  which 
is  within  5%  of  the  exact  value  for  the  10°  cone  and  for  the  15°  cone  up  to  M  =  2.5.  This 
is  considerably  better  than  the  other  methods  shown. 

(b)  Hewitt's  axisyrmnetric  body 
1 2 

Recently  the  wave  drag  of  a  particular  axisymmetric  body  has  been  calculated 
using  both  the  method  of  characteristics  and  the  exact  solution  to  the  linear  equation. 

The  radius  of  the  body  in  terms  of  the  streamwise  coordinate  x  ,  is  given  by 

R  =  Jx(l  -  x)^  0  <  x  <  1  (23) 

where  the  body  is  assumed  to  be  of  unit  length  and  has  a  cusp  at  its  downstream  end 
(ie  at  x  =  1).  The  shape  of  the  body  and  the  wave  drag  from  both  the  method  of  character¬ 
istics  (continuous  line)  and  the  exact  solution  of  the  linear  equation  (dotted  line)  are 
shown  in  Fig  10.  The  good  agreement  between  these  values  is  not  matched  by  the  super¬ 
sonic  area  rule  estimate,  as  shown  by  the  dashed  line  for  N  =  N  =100. 

•'  y  z 

We  saw  in  the  case  of  the  cone-cylinder  that  the  supersonic  area  rule  estimate  was 
close  to  the  exact  value  for  cone  angles  of  less  than  about  half  the  Mach  angle,  but  that 
large  over-estimates  of  the  drag  were  likely  to  occur  for  larger  cone  angles.  For 
M  ^  1.6  the  nose  angle  of  the  body  of  equation  (23)  is  more  than  2/3  of  the  Mach  angle, 
and  hence  a  high  wave  drag  estimate  is  not  necessarily  in  conflict  with  the  cone-cylinder 
comparison.  A  more  detailed  investigation  of  the  error  is  necessary  therefore,  to  estab¬ 
lish  whether  the  smooth  axisymmetric  bodies  of  equation  (23)  have  an  error  pattern  con¬ 
sistent  with  that  of  cone-cylinder  bodies. 

First,  consider  the  wave  drag  of  bodies  which  are  described  by  equation  (23),  but 
which  are  truncated  such  that  the  base  is  not  at  x  =  I  but  at  some  point  0  <  x^^  <  1 . 

Thus  when  x^^  is  small  the  body  resembles  a  cone-cylinder  body  and  when  x  =  1  the  body 
is  that  defined  by  equation  (23).  The  wave  drag  at  M  =  1.6  for  these  bodies  is  shown 
in  Fig  11,  both  for  SAR  and  the  method  of  characteristics,  and  it  can  be  seen  that  the 
large  differences  between  them  found  at  x^^  =  1  in  Fig  10,  occur  for  most  of  the  range  of 
x^  .  The  error  in  the  SAR  estimate,  when  expressed  as  a  percentage  of  the  method  of 
characteristics  value  is  shown  in  Fig  12.  The  variation  in  the  error  over  the  forebody 
and  the  rather  higher  value  established  over  the  rest  of  the  body  suggests  that  the 
large  nose  angle  only  partially  explains  the  drag  error. 

Cone-cylinder  results  would  suggest  that  the  error  in  the  wave  drag  should  decrease 
rapidly  as  the  body  becomes  more  slender.  In  Fig  13  the  wave  drag  for  bodies  of  unit 
length  whose  radius  is  given  by 

“  Xx(l  -  x)^ 


R 


0  <  x  <  1 


(24) 


M 


are  shown.  As  expected  the  error  decreases  as  the  body  becomes  more  slender  (X  smaller), 

but  even  for  X  <  0.35,  when  the  nose  angle  is  less  than  half  the  Mach  angle,  errors  of 

ISZ  in  the  wave  drag  estimate  using  SAR  still  occur.  This  is,  however,  significantly 
better  than  the  estimate  using  normal  cuts,  which  is  shown  in  Fig  13  to  have  an  error 
of  nearly  25%  when  X  «  0.35. 

The  lack  of  consistency  in  the  results  even  for  these  comparatively  simple  axi- 
synmetric  shapes  indicates  that  comparisons  are  needed  for  a  much  wider  range  of  shapes 
before  any  more  general  conclusions  can  be  reached.  The  encouraging  factor  in  the 
examples  shown  is  the  significant  improvement  in  the  drag  estimates  from  the  'oblique 
cuts’  method  used  here  when  compared  with  the  estimates  from  'normal  cuts'.  It  is 

hoped  that  the  provision  of  a  simple  technique  for  obtaining  oblique  cross-sectional 

areas  will  encourage  wave  drag  estimates  based  on  oblique  cuts  to  be  included  in  future 
wave  drag  comparison  studies,  enabling  the  strengths  and  weaknesses  of  this  method  to 
be  compared  with  those  of  other  methods. 

5  CONCLUSIONS 

A  new  method  of  obtaining  oblique  cross-sectional  areas  of  general  body  shapes  is 
shown  to  be  both  efficient  and  accurate.  When  used  in  conjunction  with  Eminton's  method 
for  evaluating  the  slender-body  wave  drag  integral,  analytical  slender-body  results  are 
reproduced  within  1%.  The  method  is  used  to  estimate  the  wave  drag  of  cone-cylinder 
bodies  and  smooth  axisymmetric  bodies  for  which  the  exact  wave  drag  is  known  from  the 
method  of  characteristics  with  a  'fitted'  bow  shock  wave.  For  cone-cylinders  with  nose 
semi-angles  up  to  half  Che  Mach  angle,  the  estimates  of  the  wave  drag  were  found  to  be 
within  5%  of  the  exact  values.  This  standard  of  accuracy  was  not  maintained  for  the 
smooth  axisymmetric  bodies,  however,  when  typically  over-estimates  of  15%  occurred  for 
nose  semi-angles  equal  to  half  the  Mach  angle.  Even  this  result  represents  a  significant 
improvement  in  the  accuracy  obtained  using  normal  cuts,  but  results  for  a  larger  range 
of  bodies  are  needed  before  any  general  conclusions  can  be  drawn. 
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THE  SPECIFICATION  OF  BODY  SHAPE  AND  CALCULATION  DATA  FOR  THE  WAVE  DRAG  PROGRAM 

The  program  is  written  for  easy  and  convenient  use.  The  body  is  defined  by  the 
Fortran  FUNCTION  X(Y,Z,IB).  This  body  function  will  in  general  be  supplied  by  the  user 
and  must  give  the  X  (streamwise)  value  of  the  body  surface  for  any  given  Y  and  Z 
value  (with  Y  =  0  a  plane  of  symmetry).  Often  there  will  be  more  than  one  value  of  X 
for  any  Y  and  Z  value,  that  is  one  value  of  X  for  each  time  this  line  enters  or 
leaves  the  body.  IB  defines  which  intersection  is  required,  such  that  IB  =  1  is  the 
first  (most  upstream)  body  entry  X  value  and  further  intersections  are  taken  in 
ascending  X  value.  For  a  simple  closed  body  X  will  have  significant  values  only 
when  IB  is  I  or  2.  An  annotated  listing  of  the  main  program  and  the  X  function 
used  appears  at  the  end  of  this  Appendix.  A  more  generally  applicable  X  function  for 
bodies  of  form  F(x,y,2)  =  0  is  also  listed. 

For  convenience,  the  run  parameters  are  defined  in  the  FUNCTION  DATA(N) .  These 
are  described  below  and  examples  of  typical  DATA  functions  are  included  in  the  program 
listing . 

RM  Mach  number  (>1 ) 

ZO  2  coordinate  of  the  X  line  about  which  the  wave  drag  is  calculated 

(value  usually  zero) 

NBODY  A  reference  body  number  for  users'  convenience 

NOPDEV  Output  channel  for  results 

NOPCTL  Controls  output  of  results 

Value  1  Run  parameters  and  wave  drag  only 

2  Wave  drag  for  each  ?  value  given  also 

3  Areas  for  each  9  and  DX  given  also 

4  Body  projection  and  XO  limits  given  also 

NX  Number  of  oblique  cuts  for  each  value  of  9  .  (NX  <  100) 

NY  Number  of  points  across  body  in  Y  direction 

(NY  <  100  but  can  be  increased  by  altering  DIMENSION  statement  in 
main  program) 

NZ  Number  of  points  across  body  in  Z  direction 

(NZ  <  100  but  can  be  increased  by  altering  DIMENSION  statement  in 

main  program) 

NTHETA  Number  of  6  values  from  0  to  it  at  which  wave  drag  is  calculated 

(NTHETA  <  100  but  can  be  increased  by  altering  DIMENSION  statement 
in  main  program) 

NB  Maximum  number  of  times  any  X  line  enters  or  leaves  the  body 

(if  NB  >  2  alter  DIMENSION  statement  in  main  program) 


(JUU  ■ 
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Appendix 


MAIN  PROGRAM. 


PROGRAM  UAVEDRAG. 

C  CALCULATES  THE  UAUEDRAG  OF  AN  ANALYTICALLY  DEFINED  BODY  X=X(YfZ) 

C  USING  SUPERSONIC  AREA  RULE  WITH  ALL  OBLIQUE  CUTS. 

C  ARRAYS  XBM(NY.N2fNB) .AREA(NX).WDRAG<NTHETA> .ZMINV(NY) »ZMAXV<NY) .DZ(NY) 
C  MAX  VALUES  OF  NY,NZ»NB.NTHETA  SET  BY  DIMENSION  »  COMMON  FOLLOWING. 
DIMENSION  XBM< 100. 100.2) .AREA( 100) .WDRAG( 100) .DZ( 100) 

COMMON/A/  ZMINU(100).ZMAXV(100) 

C  ZMINV  S  ZMAXV  ARE  ONLY  IN  COMMON  FOR  THE  SIZE  CHECK  OF  ZMINV  NY  BELOW 

C 

c  GET  RUN  DATA. 

RM=DATA(1) 

Z0=DATA(2) 

N0FDEM=DATA(3)+0- 1 
N0PCTL=DATA(4 '+0. 1 
NBODY-DATA' 5>+0. 1 
NX.=  DATA(6)+0 . 1 
NY~DArA(7)+0. 1 
NZ=DATA<8)+0. 1 
NTHElA=DATA<<?i+0. 1 
NB=DATA<10)+0.1 

0 

r.  CHECK  NY  SIZE  AGAINST  ZMINV  SIZE  ALLOTTED  IN  COMMON  STATEMENT- 
DO  1990  N=1.NY 
1990  ZMINV<N)=0 

ZMAXV< 1 >=99999 
DO  2000  N=2.NY 

2000  IF<ZMINV(N) . GT . 99998 . 99 . AND . ZMINV( N ) .LT. 99999. 01 /GOTO  2010 
G0T02030 

2010  WRITE(N0PDEV.2020> 

2020  FORMAT<"  NY  BIGGER  THAN  ARRAY  SPACE  GIVEN  IN  DIMENSION  STATEMENT'’ 
WR I TE  <  NOPDEV . 1 10 ) NX . NY . NZ . NTHETA . NB 
GOTO  1100 
2030  CONTINUE 
C 

f  DESCRIBE  THE  BODY  PROJECTION  SHAPE  <0N  A  PLANE  X=CONST.> 

C  INITIAL  GUESS 
YMAX=0.1 
ZMIN=-0. 1 
ZMAX=0.1 
C 

C  FIRST  EXPAND  YMAX.ZMIN  S  ZMAX  TO  CONTAIN  PROJECItOH. 

1  YMAXO=YMAX 
ZMINO=ZMIN 
ZMAXO=ZMAX 
DO  2  IY=1.NY 
Y=<IY-0.5)*YMAX/NY 

IF ( X ( Y . ZMAX . 1 ) . LT . 9000 ) ZMAX=1 . 2#ZMAX-0 . 2*ZM IN 
IF<X<Y.ZMIN.l) .LT.9000)ZMIN=1.2*ZMIN-0.2»ZMAX0 
IF(ZMAX-ZMAXO+ZMINO-ZMIN.GT.  (ZMAX-ZMIN)  ■'20)G0T0  I 

2  CONTINUE 

DO  3  IZ=1.NZ 

Z=ZMIN+(IZ-0.5)#(ZMAX-ZMIN)/NZ 
I F  <  X  <  YMAX . Z . 1 ) . LT . 9000 ) YMAX= 1 . 2»YMAX 
IF(YMAX-YMAX0.GT.YMAX/20)G0T0  1 

3  CONTINUE 
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C 

NOW  REUUCt:  iflA)- »ZMiN»2MPtX  TO  HUG  PROJECTION 
C  NEW  ZMIN  ii  ZMAX  VECTORS  IN  ZNINV^NT)  f  ZMAXV(NY) 

nzo:-  ■>  ZhAX-ZhiN )  'Nz/s 

Y  =  t'MAX 
YNAX-r 

DO  •’  Ti'-1,NY 
ZMAX  "ZiiAXO 
ZMIN-XZMINO 
I  :  Y  -NYM  -lY 

Y  =  •'  1 1 Y  -0 . 5  )  #  YMAX  /  N  Y 
Z--ZMAX+DZ0/2 

5  : -Z-DZO 
IF<Z.L.E,ZMIN)GOTO  4 
IP(X<YfZfl).GT.9000)G0T0  5 
ZMAXV( IIy)=AMIN1 (ZMAXfZ+DZO) 

Z=ZMIN-riZ0/2 

6  Z“Z+jzo 

IF  (X<  Y»Z»1)  .GT.9000)G3T0 
ZhINVI IIY)=AhAXl(ZMIN»Z-DZO) 

7  CONTINUE 

r; 

c;  set  up  IiY,DZ!.Ii  THETA  S  A  GUESS  AT  XSTARU  XF  INIS . 

XS--9000 

Xf-::::-9000 

IF<  ABSIXIOfO,  1  )  -  .l.T,9000)XS=X(0,0»l  ) 

IF"<AE(S<X(0»0,NB)  )  ,L.T.9000)XF-=-X(0f0.Nro 
riY-:YMAX/  NY 
DO  8  I Y^^  It  NY 

8  DZ  ■:  I  Y  (  ZMAXV  ( I Y  )  -ZMINV<  I Y  )  )  /NZ 
0THETA-=:3,  1.4159/  <NTHETA-1  > 

r 

C  SET  U(-  ./  value  AF'PAYS  (FOP  PROJECTED  SHAPE.) 

DO  10  :i:i^=-l.ND 
DO  10  IY^:=lrNY 
DO  10  IZ^=1tNZ 
XDM( lY  > IZt IB) -20000 
Y^^aY-0.5)*DY 

Z~ZM I NV  (  I  '^ )  f  (  1 2-0 . 5 )  *DZ  (  I Y  ) 

IIB=^  IB 

XVALUE^=X(Y,Zt1IB> 

XBM(IYTlZrIB)=XVALUE 

IF( ABS<XVALUE) .GT, 9000) GOTO  10 

C  FIND  UPSTF'EAM(XSTART)  &  DOWNSTREAM(XFINIS)  POINTS  ON  REF  AXIS! Y=OtZ=ZO) 
C  FOF(  WHICH  D(  OBLIQUE  AREA)/DX=0. 

RBETA--:SQRT<  (  Y)«  YK  Z-ZO  )  *  ( Z-ZO  )  J’KRMjKRM-I  )  ) 

XS-AMINl ( XS  T  XVAl UE~RBETA ) 

XF::::AMAX1  <  XF  T  XVAl.UE  +  RBETA  > 

10  CONTINUE 

XSTART::^XS-(XF-XB)  NX 
XFINIS==XF  f  (  XF  -XS  :)/NX 
DX:::(  XF  INIS  -XSTART)  -'(NX-i  .  > 


( 

i 
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i^pendix 


URITE  RUN  DATA. 

WRITE  <  NOPDEUf 100 ) RM  f  ZO  »NBODY / NOPDEV » NOPCTL 
100  FORMAT(/'  DATA'>/»'  MACH  NUMBER<RM)=' rF6.4» ' f  AXIS  Z(ZO)=' fF5.3f/f 
8'  NB0DY='Fl2f'  NOPDEU=' fI2f ' F  NOPCTL= ' f 12) 

URITE  <  NOPDEV  f 1 10 ) NX  f  NY  f  NZ  f  NTHETA  f  NB 
110  FORMAT<'  NX='fI3f'f  NY='fI3f'f  NZ='fI3f'f  NTHETA='  fI3f'f  NB='fI3) 
IF(N0PCTL.LT.4)G0T0  20 
WRITE ( NOPDEO  f 1 20 ) XSTART  f  XFINI S  f  YMAX 
120  FORMAT</'  XSTART=' fF8.3f ' f  XFINIS=' fF8.3f ' f  YMAX='fF8.3) 
WRITE<N0PDEUf130) 

130  FORMAT('  ZMINVfZMAXV' ) 

WRITE<N0PDEVf140)<ZMINV< J) F J=1fNY> 

WRITE (NOPDEOf 140 ) (ZMAXV(J) fJ=1f NY) 

140  F0RMAT(1Xf10F7.4) 

20  CONTINUE 


C  EVERYTHING  IS  NOW  SET  UP  SO . 

C  FIND  OBLIQUE  AREAS  8  WAVEDRAG  FOR  EACH  THETA. 

C 

i;  MAIN  LOOP  (THETA) 

DO  60  ITHETA=1fNTHETA 
THETA=DTHETA»( ITHETA-1 ) 

BS=SaRT<RM»RM-l )*SIN<THETA) 

BC==SQRT  ( RMiKRM- 1 )  »COS  <  THETA ) 

C 

C  2ND  LOOP  <X> 

DO  50  IX=1fNX 
XO~XSTART+<  IX-1  )ilcDX 
SUMKUT=0 

C 

I-  3RD  LOOP  (Y) 

NYY^:=1-NY 

DO  40  IY=NYYfNY 

IIY=<IABS<2*IY-l)+l)/2 

Y:=(IY-0,5)*DY 

PARTYZ=XO+BS»Y+BC*<ZMINV<IIY)-DZ( IIY)/2-Z0) 
PARTZ=BC*DZ<IIY) 

C 

C  4TH  LOOP  (B) 

DO  40  IB=1fNB 

DZIIY=(2*M0D<IBf2)-1  ))HDZ(IIY) 

KUT0LD=0 

CUTOLD=20000 

C 

C  5TH  LOOP  (Z) 

DO  40  IZ=1fNZ 
KUT=0 

XCUT=PARTYZ+IZ»PARTZ-XBM< IIYf IZf IB) 
IF(XCUT.GT.0)KUT=1 

XCUTS=AMIN1< 1 ,0fAMAX1(0.0fXCUT/DY/2.0+0.5) ) 
SUMKUT=SUMKUT+XCUTSi|iDZIIY 
IF<KUT,EQ.KUTOLD)GOTO  30 


on  n  no 
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INTERPOLATE  IN  Y  PLANE  IF  SURFACE  POINT. 

XCUTY=XCUT 
XCUTOY=CUTOLD 

DENOM=ABS  ( XCIIT )  +ABS  <  CUTOLD ) 

IF (DENOM.lt. 10000) GOTO  25 
EXTRAPOLATE  TO  EDGE  OF  NON-POINTED  BODY. 

IF<ABS(CUTOLD) . GT. 10000. AND. ABS(XCUT) .GT. 10000)G0T0  25 
PART=PARTYZ+IZ»PARTZ-2*XBH( IIY» IZ» IB) 

IF (ABS( CUTOLD) .GT. 10000)XCUT0Y=PART+XBM( IIY» IZ+1 r IB) 

IF ( ABS ( XCUT ) . GT .10000 ) XCUTY=PART+XBM (IIY»IZ-1,IB) 

DENOM=ABS  <  XCU  TY ) +ABS  <  XCUTOY ) 

25  SUMKUT=SUMKUTf ( XCUTY+XCUTOY ) /DEN0M/2*DZI I Y 
30  KUTOLD^KUT 
CUTaLD==XCUT 
40  CONTINUE 

NOU  CALC.  AREA  S  FROM  AREAS  GET  WAVE  DRAG  FOR  THIS  THETA 
AREA< IX)=SUMKUT*DY 
50  CONTINUE 

CALL  WAUEDR<DX.AREA>NXrUDRAG(ITHETA) ) 

THETAD=THETA»180/3. 14159 

IF^MODdTHETA-l.l )  .  NE.  0.  AND.  I  THETA.  NE.NTHETA)  GOTO  60 
IF (NOPCTL.lt. 2) GOTO  60 
UR I TE ( NOPDEU » 200 ) THETAD  »  UDRAG ( I THETA ) 

200  FORMAT(/'  THETA=' »F6.2. ' ,  WAVEDRAG=' .F12.6) 

IF (NOPCTL.lt. 3) GOTO  60 
URITE(N0PDEV»210) 

210  FORMAT('  NX  AREAS') 

WRITE(N0PDEU»220) ( AREA( J ) f J=1 » NX ) 

220  F0RMAT(1X»10F7.4) 

60  CONTINUE 

C  AVERAGE  UAVEDRAG  OVER  THETA  (AWDRAG) 

UDRAG(l) -UDRAG (l)/2 

UDRAG ( NTHETA ) =UDRAG ( NTHETA ) /2 

AUDRAG=0 

DO  70  ITHETA=1, NTHETA 
70  AUDRAG-AUDRAG+UDRAG ( ITHETA ) / ( NTHETA- 1 ) 

UR I TE ( NOPDEV . 300 ) AUDRAG 

300  FORMATv//'  BODY  UAVEDRAG=' fF13.7F '  (SQ.  LENGTH  UNITS  OF  X»Y»Z)') 
1100  CONTINUE 
STOP 
END 


in 
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SUBROUT I NE  UAUEDR ( DX  f  AREA  »  NX » WDRA6 ) 

C  NX  AREAS  ARE  SUPPLIED  IN  AREA>SPACED  AT  INTERVALS  DX 
C  UDRAQ  IS  PUT  EQUAL  TO  THE  UAVE  DRAG  (D/Q)  USING  EHINTONS  METHOD 
0 1 MENS I ON  AREA ( NX ) f  C ( 1 00 ) f  P  < 1 00  » 1 00  >  f  WORK 1 ( 1 00 ) f  U0RK2 ( 1 00 ) 

NX2»NX-2 
DO  10  1X^=1  fNX2 
X® IX/ ( NX— 1 • ) 

Us(2*ATAN(SQRT(X/<l-X) ) )-2*< 1-2#X>»SQRT ( X-X*X ) >/3 . 14159 
C <  I X > » AREA (  I X-f  1 ) -  AREA <!)-< AREA < NX ) -  AREA (1))*U 
DO  10  IY==lfNX2 
Y=IY/(NX-1. ) 

XY1=X+Y-2*X*Y 

XY2=2*SQRT<X*Y*<1-X)*<1-Y) ) 

P<IX»IY)=XY1*XY2 

IF(IY.NE.IX)P<IX»IY)=-<X-Y)*»2»AL0G( (XY1+XY2)/(XY1-XY2) )/2+XYl»XY2 
10  CONTINUE 

C  INVERT  THE  ARRAY  Pa00,100)  ONTO  ITSELF 
CALL  I NVERT  <  P  f  NX2 . 1 00  f  WORK 1 »  W0RK2 ) 

SUM=0 

DO  20  IX=1«NX2 
DO  20  IY=1»NX2 

20  SUM=SUM+C(IX)*C<IY)#P(IX.IY) 

UDRAG=<4/3.1416»(AREA<NX)-AREA<1) >»*2+3.1416#SUM>/DX/DX/(NX-l)**2 
RETURN 
END 
C 
C 

c 

C; 

SUBROUTINE  INVERT< Af M> lAr IND»C) 

C  INVERTS  AN  M»M  MATRIX  ‘A*  ONTO  ITSELFfWHERE  AI  IS  THE  DECLARED 
C  ARRAY  SIZE  AND  IND<M)  S  C<M)  ARE  WORKING  ARRAYS. 

C 

C  LIBRARY  SUBROUTINE  INVERT  CODE  NOT  INCLUDED, 


v,n 

ijt 
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C  BODr  t  DATA  FUNCTIONS  FOR  SEARS-HAACK  BODY. 

C 

C 

FUNCTION  X<Y»Z»IB) 

C  SEARS-HAACK  BODY  OF  GIUEN  LENGTH  (XL)  t  D/Q. 

C  WHEN  M=1.001  THIS  TEST  BODY  SHOULD  GIVE  D/Q  FROM  SAR  THE  SAME  AS 
C  THAT  ENTERED  BELOW  (FOR  ANY  AOB  t  LENGTH  ) 

C  NB=2  FOR  THIS  BODY. 

C  RMAX=(2«L«L«D0Q/9/RI««3>t«(l/4>  PG283  NIELSENf MISSILE  AERODYNAMICS 
C  AOB=RATIO  A/B  FOR  ELLIPTIC  CROSS-SECTION  OPTION. 

X=20000 

A0B=1.0 

XL=1 

D0Q=0.08 
PI=3. 14159 

TEMP=9*PI#PI)|tPI* ( ( Y*Y/AOB+Z#Z*AOB )  **2 ) /2/XL/XL/DOQ 
IF(TEMP.GT.1.0R.IB.GT.2)G0T0  10 
X=XL/2*(  1 -SORT  <  l-TEMP)lc*(  1/3. ) ) ) 

IF(IB.EQ.2)X=XL-X 
10  RETURN 
END 
C 
C 


FUNCTION  DATA(N> 

RM=1.0001 

Z0=0 

N0PDEV=5 

N0PCTL=4 

NB0DY=1 

NX=25 

NY=50 

NZ=50 

NTHETA=3 

NB=2 

IF(N.Ea.l)DATA=RM 

IF(N.EQ.2)DATA=Z0 

I F ( N . EQ . 3 ) DATA=NOPDEU 

I F ( N . EQ . 4 ) DATA=NOPCTL 

IF ( N . EQ . 5 ) DATA=NBODY 

IF(N.EQ.6)DATA=NX 

IF(N.EQ.7)DATA=NY 

IF(N.Ea.8)DATA=NZ 

IF ( N . EQ , 9 ) DATA=NTHETA 

IF(N.EQ,10)DATA=NB 

END 
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C  BODY  t  DATA  FUNCTIONS  FOR  SINGLE  BULGE  BODY. 

C 

C 

FUNCTION  X<Y»Z»IB) 

C  SINGLE  BULGE  BODY  (NB^Z  FOR  THIS  BODY) 

X»20000 

IF<ABS<Z>.GT,1.5)Q0T0  999 
IF<ABS(Y).GE.2.5)G0T0  999 
IFdB.GT.DQOTO  10 
IF<ABS<Y) .LE.1.S)X=-20000 
IF<X.LT.-9999)G0T0  999 
IF<Z.GE.0)X=A/3.14159*AC0S<4-2*Y) 

G0T0999 

10  IF<ABS(Y).LE.0.5.0R.ABS<Y>.6E.1.5.0R.Z.GE.0)G0T0  999 

X=A/3. 14159*AC0S<2*Y-2) 

999  RETURN 

END 
C 
C 
C 

c 

FUNCTION  DATA<N) 

RM=1.4 

Z0=0 

N0PDEV=5 

N0PCTL=4 

NB0DY=2 

NX*31 

NY*100 

NZ=100 

NTHETA=19 

NB=2 

IF(N.EQ.1)DATA=RM 

IF(N.EQ.2)DATA=Z0 

IF  <  N . EQ . 3 ) DATA=NOPDEV 

IF  <  N . EQ . 4 ) DATA=NOPCTL 

IF  <  N . EO . 5 ) DATA=NBODY 

IF<N.Ea.6)DATA=NX 

IF<N.Ea.7)DATA=NY 

IF<N.EQ.8)nATA=NZ 

I F  <  N . EQ . 9 ) DATA=NTHETA 

IF(N.Ea.lO)DATA=NB 

END 


fj  n  o  n  n  n  n  n  n  o  on 


Appendix 


21 


BODY  FUNCTION  FOR  ANY  BODY  F(X»Y»Z)=0. 


FUNCTION  X(Y»Zi-IB) 

X  VALUE  OF  ANY  FUNCTION  F<X.Y,Z>=0 
IB=1  SMALLEST  X  INTERSECTION  VALUE 


GIVEN  Y  t  Z. 

<>XN0SE).IB-2  NEXT  INTERSECTION 


USER  DEFINED  FUNCTION  F 
<BODY  OF  FIG  10  HERE'. 

F<X,y,Z)=<Y*'i'  +Z#Z  >  '0 . 25*X**2*  (  1  -X  )  **4 


USER  SET  RUN  CONSTANTS 
0  XERR0R=0, 000001 
XN0SE=0 
XBASE=1 

PROGRAM (NO  USER  CHANGE  NEEDED). 

C;  INITIALISE 
IMAX-100 
XN-XNOSE 
XB-XBASE 
1  BC=--^0 

FN=F  <0,0 f 1 0000 ,0rl 0000 . 0 ) 

C  FIND  X  FOR  IB<X  IN  X  TO  X+DX  HERE) 

C  F  CHANGES  SIGN  AS  X  CROSSES  SURFACE. 

100  DO  200  I=1.IMAX 

DX=(XB-XN)/<IMAX-1 ) 

X=XN+(I-1)*DX 
F0LD®FN 
C  GET  F  VALUE 

FN =r  r  X . Y » Z ) 

110  IF(FN»F0LD.LE.0)IBC=IBC+1 
IF<IBC.EQ.IB)G0T0  300 
200  CONTINUE 

C  NO  VALUE  FOR  THESE  Y»Z»IB  SO  SET  X=DEFAULT . 
X®20000 
GOTO  400 

300  I F  <  X . LT . XN+0 . 00000000 1 ) X=-20000 

IF<X.LT.-19999)G0T0  400 
C  CRAMP  X-LIMITS  XNSXB  ONTO  X 
XN=X-DX 
XB=X 

I6C=IBC-1 

FN=F<0,0» 10000. Of  10000.0) 
IF<DX.GT.XERROR*2)GOTO  100 
X=(XN+XB)/2.0 
IF  <  X . GT . XBASE ) X=20000 
IF  <  X . LT . XNOSE ) X=-20000 
400  RETURN 
END 
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LIST  OF  SYMBOLS 

cross “Sectional  area  of  body  element 
wave  drag 

source  distribution  function 

body  shape  function  (equations  (17)  to  (22)) 

body  length  (or  streamwise  length  for  which  S  varies) 

Mach  nianber 

number  of  cro8S“8ectional  areas 

number  of  body  elements  in  y  and  z  directions 
kinetic  pressure 
body  radius 

body  cross-sectional  area  (may  be  oblique) 
free-stream  velocity 
body  volume  (equation  (II)) 
streamwise  coordinate 

body  element  intersections  with  body  surface 
spanwise  coordinates 
coordinates  of  elemental  body 
(M^  -  1)^ 

interval  between  cross-sections 
elemental  body  spacing 
cylindrical  polar  angle 
free-stream  density 

coordinates  parallel  and  normal  to  base  contour  (equation  (D) 

thickness  parameters 

slender-body  potential  (equation  (D) 

angle  between  the  normal  to  the  cutting  plane  and  the  x  axis  (cos 


'  ( I  /M)  ) 


23 


So.  Author 

1  G.N.  Ward 


R.T.  Whitcomb 


3  J.S.  Nielsen 

4  C . S .  Leyman 
T.  Markham 

5  £.  Efflinton 

6  W.R.  Sears 

7  W .  Haack 

8  G.  Norman 

9  J.L.  Sims 


10  M.J.  Lighthill 

11  J.B.  Broderick 

12  B.L.  Hewitt 


REFERENCES 

Title,  etc 

Linearised  theory  of  steady  high-speed  flow. 

Cambridge  Monographs  on  Mechanics  and  Applied  Mathematics, 
Cambridge  University  Press  (1955) 

A  study  of  the  zero-lift  drag-rise  characteristics  of  wing-body 
combinations  near  the  speed  of  sound. 

NACA  TR  1273  (1956) 

Missile  aerodynamics. 

McGraw  Hill,  New  York,  pp  283-300  (1960) 

Prediction  methods  for  aircraft  aerodynamic  characteristics. 
AGARD  LS-67,  paper  5,  May  1974 

On  the  minimisation  and  numerical  evaluations  of  wave  drag. 

RAE  Report  Aero  2564  (ARC  R  &  M  3341)  (1958) 

On  projectiles  of  minimum  wave  drag. 

Quart.  Applied  Machs,  Vol.l4.  No. 4  (1947) 

Geschossenformen  kleinsten  Wei lenwider stands,  Ber. 
Lilienthal-Ges.  Luftfahrt,  Vol.139,  pp  14-37  (1948) 

Some  calculations  of  the  wave  drag  of  bluff  asymmetries. 

British  Aerospace,  Kingston  -  Brough  Div,  Note  YAD  3342, 
November  1 97  8 

Tables  of  supersonic  flow  about  right  circular  cones  at  zero 
angle  of  attack. 

NASA  SP  3004  (1964) 

Higher  approximations  in  aerodynamic  theory. 

Princeton  University  Press,  p  126  (1960) 

Supersonic  flow  past  a  semi-infinite  cone. 

Quart. J.  Mech.  and  Appl.  Math.  2,  pp  121-128  (1949) 

The  calculation  of  wave  drag  in  steady  supersonic  flow. 

British  Aerospace,  Warton  Division,  Report  AE/A/651,  July  1980 


Fig  4  Error  in  Mtimating  tho  wav*  drag  of  a  Saart-Haack  body 
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Fig  10  Drag  oowfficitnt  of  a  ci<md  axifynwimric  body 
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Fig  11  W«v*  drag  of  a  smooth  axiiymnwtrie  body  truncattd  at  x  -  X|_ 
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Fig  13  Wave  drag  of  a  HiMOth  Mitymmetric  body  for  •  range  of  thidcnMi 


